指数関数のパラメータの推定
a0,をシフト
\(\Large \displaystyle y_i = a_0 \ Exp (- a_1 x_i) + a_2 \)
a0 | 9.1640 | 9.6640 | 9.9640 | 10.0640 | 10.1640 | 10.2640 | 10.3640 | 10.6640 | 11.1640 | |||
δ | -1.0 | -0.5 | -0.2 | -0.1 | 0.0 | 0.1 | 0.2 | 0.5 | 1.0 | |||
a1 | 0.5321 | 0.5134 | 0.5021 | 0.4983 | 0.4946 | 0.4908 | 0.4870 | 0.4755 | 0.4559 | |||
a2 | 1.3123 | 1.1040 | 0.9735 | 0.9290 | 0.8839 | 0.8383 | 0.7921 | 0.6500 | 0.3993 | |||
i | x | y | \( \hat{y} \) | |||||||||
1 | 0 | 10 | 10.4763 | 10.7680 | 10.9375 | 10.9930 | 11.0479 | 11.1023 | 11.1562 | 11.3140 | 11.5633 | |
2 | 2 | 4 | 4.4739 | 4.5654 | 4.6236 | 4.6436 | 4.6639 | 4.6846 | 4.7055 | 4.7702 | 4.8850 | |
3 | 3 | 2 | 3.1694 | 3.1756 | 3.1828 | 3.1858 | 3.1891 | 3.1928 | 3.1968 | 3.2111 | 3.2427 | |
4 | 4 | 1 | 2.4031 | 2.3438 | 2.3107 | 2.3001 | 2.2897 | 2.2796 | 2.2698 | 2.2419 | 2.2017 | |
5 | 6 | 0.5 | 1.6886 | 1.5481 | 1.4633 | 1.4350 | 1.4067 | 1.3784 | 1.3501 | 1.2651 | 1.1235 | |
6 | 9 | 0.1 | 1.3886 | 1.1992 | 1.0821 | 1.0425 | 1.0025 | 0.9622 | 0.9216 | 0.7977 | 0.5837 | |
S (\(y_i - \hat{y} \)の平方和) | 0.8610 | 0.4039 | 0.2771 | 0.2591 | 0.2531 |
0.2591 | 0.2769 | 0.4010 | 0.8384 | |||
dS (Seとの差分) | 0.6078 | 0.1507 | 0.0240 | 0.0060 | 0 | 0.0060 | 0.0238 | 0.1479 | 0.5853 | |||
・残差平方和
推定値からの残差
\(\Large \displaystyle Se = \sum_{i=1}^{n} \left[ y_i -\hat{a_0} \ Exp(- \hat{a_1} \ x_i) - \hat{a_2} \right]^2 \)
a0をシフトさせたときの,推定値からの残差
\(\Large \displaystyle Se = \sum_{i=1}^{n} \left[ y_i -a_0 \ Exp(- \hat{a_1} \ x_i) -\hat{a_2} \right]^2 \)
であり,a0を,δ,だけシフトさせて,固定し,その際のa1, a2の推定値をソルバーで推定しました.
dS,を見ていただけるとわかるように,推定値,Seが一番小さく,ほぼ左右対称に増加していることがわかります.
グラフ化すると,
のように,二乗+定数できれいに近似できます.
二次曲線の近似においてもきれいに近似でき,
\(\Large \displaystyle y = 0.2532 + 0.5965 \ \delta^2 \)
ここで,分散値は,
・分散
\(\Large \displaystyle Ve = \frac{1}{n-3} \sum_{i=1}^{n} \left[ y_i -\hat{a_0} \ Exp(- \hat{a_1} \ x_i) \right]^2 = \frac{Se}{n-3} = \frac{0.2532}{6-3} = 0.08438 \)
であり(a0,a1,a2の二つのパラメータが3つあるので,自由度は,n-3),
\(\Large \displaystyle 0.5965 \ \delta^2 = 0.08438 \)
となるδがSEとなるので,
\(\Large \displaystyle \delta^2 = \frac{ 0.08438}{0.5965} =1.4145 \)
\(\Large \displaystyle SE_{a_0} = \sqrt{\delta^2} =\color{red}{0.3761} \)
と推定できます.
・Rによる推定
Rでの近似を行ってみると,
プログラムは,
xx <- c(0,2,3,4,6,9)
yy <- c(11,5,3,2,1.5,1.1)
plot(xx,yy)
fm<-nls(yy~a0*exp(-a1*xx)+a2,start=c(a0=10,a1=0.5,a2=1),trace=TRUE)
summary(fm)
で,結果は,
Parameters: | |||||
Estimate | Std. Error | t value | Pr(>|t|) | ||
a0 | 10.16402 | 0.37708 | 26.954 | 0.000112 | *** |
a1 | 0.49456 | 0.04408 | 11.219 | 0.001518 | ** |
a2 | 0.88393 | 0.26931 | 3.282 | 0.046348 | * |
となり,Kyplotにおいても,
推定値 | 標準誤差(SE) | |
A1 | 10.16401 | 0.377081 |
A2 | 0.494562 | 0.044083 |
A3 | 0.883926 | 0.269308 |
と同じ結果となり,今回の推定値と,ほぼ一致,します.
微妙に異なるのが気になりますが....
「統計解析の初歩」,の「1.2 非線形最小2乗法の基本的な考え方」には,
δのずらした値0.5 を変えると結果は異なり、近似標準誤差の精度が変わる
統計パッケージにより、近似標準誤差の値は幾分異なる
とあります.ここで,”0.5”,がどこから出てきたかはわかりません.そもそも横軸(x軸)の範囲に依存しちゃいますし..
そこで,σをとる値(±での平均値)でどう推定値が変わるかを調べてみました.その結果が,
とどのδでもRなどの推定値より下回っていました....謎です...
次に,指数関数+baseの肩のパラメータ,a1について,確認してましょう.